%% demo with MLR-GCK with AVIRIS Indiana Pines dataset 
%% for the generalized compposite kernel paper intoduced in the following paper
%
%%
%  [1] J. Li, P. R. Marpu, A. Plaza, J. M. Bioucas-Dias and J. A.Benediktsson. 
%  Generalized Composite Kernel Framework for Hyperspectral Image Classification.
%  IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 9, pp. 4816-4829,2013.
%
%%
% If you use this demo, please be aware of the parameters involved in the
% building of the kernel functions and in the LORSAL algorithm, which might
% affect the final performance the GCK framework. 
%
%%
% It is also recommended to use other spatial features for evaluation
%
%%
% For any question and suggestion, it is greatly appreciated to contact the
% authors: jun@lx.it.pt
%

close all
clear all
clc

load EMP_Indina_profiles325 % EMP data
sz = size(EMP_KPCA_Indian);
EMP_KPCA_Indian = reshape(EMP_KPCA_Indian, [prod(sz(1:2)) sz(3)]);   
EMP_KPCA_Indian = EMP_KPCA_Indian';

load imgreal; % original data
sz = size(img);
img = reshape(img, [prod(sz(1:2)) sz(3)]);   
img = img';
img([104:108 150:163 220],:) =[]; % remove noise bands


% load ground truth data
load AVIRIS_Indiana_16class
trainall = trainall';


indexes=train_test_random_new(trainall(2,:),10,160);
train1 = trainall(:,indexes);  % training sample
test1=trainall;                % remaining samples are used for test
test1(:,indexes)=[];
        
no_classes=16;
    
train_spe = img(:,train1(1,:));
[d,n] =size(train_spe);
nxo = sum(train_spe.^2);
[X,Y] = meshgrid(nxo);
disto=X+Y-2*train_spe'*train_spe;
scaleo = mean(disto(:));
clear X Y nx
sigmao = 1.5;        
% build design matrix (kernel)
Kspe=exp(-disto/2/scaleo/sigmao^2);        

train_spatial = EMP_KPCA_Indian(:,train1(1,:));
[d,n] =size(train_spatial);
nxi = sum(train_spatial.^2);
[X,Y] = meshgrid(nxi);
disti=X+Y-2*train_spatial'*train_spatial;
scalei = mean(disti(:));
clear X Y nx
sigmai = 2.0;        
% build design matrix (kernel)
Kspa=exp(-disti/2/scalei/sigmai^2);
        
K = [ones(1,n); Kspe;Kspa]; 
        
tic
lambda = 0.00001;
[w,L] = LORSAL(K,train1(2,:),lambda,400);

toc
%%     GEK
p=splitimage_MLR_GCK(img,EMP_KPCA_Indian,train_spe,train_spatial,w,scaleo,scalei,sigmao,sigmai);
[maxp,class] = max(p);

[overall,kappa,average,classindividual] = calcError(test1(2,:)-1,class(test1(1,:))-1,[1:no_classes]);
        
overall
